Thursday, September 21, 2022
Today
- Making maps
Making maps using functions from the {tmap} package
The {tmap} package has functions for creating thematic maps. The syntax is like the syntax of the functions in {ggplot2}. The functions work with a variety of spatial data.
Consider the simple feature data frame called World from the {tmap} package.
library(tmap)
data("World")
str(World)## Classes 'sf' and 'data.frame': 177 obs. of 16 variables:
## $ iso_a3 : Factor w/ 177 levels "AFG","AGO","ALB",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ name : Factor w/ 177 levels "Afghanistan",..: 1 4 2 166 6 7 5 56 8 9 ...
## $ sovereignt : Factor w/ 171 levels "Afghanistan",..: 1 4 2 159 6 7 5 52 8 9 ...
## $ continent : Factor w/ 8 levels "Africa","Antarctica",..: 3 1 4 3 8 3 2 7 6 4 ...
## $ area : Units: [km^2] num 652860 1246700 27400 71252 2736690 ...
## $ pop_est : num 28400000 12799293 3639453 4798491 40913584 ...
## $ pop_est_dens: num 43.5 10.3 132.8 67.3 15 ...
## $ economy : Factor w/ 7 levels "1. Developed region: G7",..: 7 7 6 6 5 6 6 6 2 2 ...
## $ income_grp : Factor w/ 5 levels "1. High income: OECD",..: 5 3 4 2 3 4 2 2 1 1 ...
## $ gdp_cap_est : num 784 8618 5993 38408 14027 ...
## $ life_exp : num 59.7 NA 77.3 NA 75.9 ...
## $ well_being : num 3.8 NA 5.5 NA 6.5 4.3 NA NA 7.2 7.4 ...
## $ footprint : num 0.79 NA 2.21 NA 3.14 2.23 NA NA 9.31 6.06 ...
## $ inequality : num 0.427 NA 0.165 NA 0.164 ...
## $ HPI : num 20.2 NA 36.8 NA 35.2 ...
## $ geometry :sfc_MULTIPOLYGON of length 177; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:69, 1:2] 61.2 62.2 63 63.2 64 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:15] "iso_a3" "name" "sovereignt" "continent" ...
The spatial data frame contains socioeconomic indicators from 177 countries around the world. Each row is one country’s indicators.
You make a map by first specifying the spatial data frame using the tm_shape() function and then you add a layer consistent with the geometry.
For example, if you want a map showing the index of happiness (column name HPI) by country, use the tm_shape() function to identify the spatial data frame World then add a fill layer with the tm_polygons() function.
The fill is specified by the argument col = indicating the specific column from the data frame. Here use HPI.
tm_shape(shp = World) +
tm_polygons(col = "HPI")
The tm_polygons() function with the argument col = colors the countries based on the values in the column HPI of the World data frame.
Map layers are added with the + operator.
Caution: the column in the data frame World must be specified using quotes "HPI". This is different from the functions in the {ggplot2} package.
To show two thematic maps together each with a different variable, specify col = c("HPI", "well_being")
The tm_polygons() function splits the values in the specified column into meaningful groups (here 8) and countries with missing values (NA) values are colored gray.
More (or fewer) intervals can be specified with the n = argument, but the cutoff values are chosen at appropriate places.
Example: Mapping tornadoes
Consider the tornado data from the U.S. Storm Prediction Center (SPC). It is downloaded as a shapefile in the directory data/1950-2018-torn-aspath.
A shapefile is imported with the sf::st_read() function from the {sf} package.
Tornadoes.sf <- sf::st_read(dsn = "data/1950-2018-torn-aspath")## Reading layer `1950-2018-torn-aspath' from data source
## `/Users/jameselsner/Desktop/ClassNotes/QG-2022/data/1950-2018-torn-aspath'
## using driver `ESRI Shapefile'
## Simple feature collection with 63645 features and 22 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -163.53 ymin: 18.13 xmax: -64.9 ymax: 61.02
## Geodetic CRS: WGS 84
The assigned file is a simple feature data frame with 63645 features (observations) and 23 fields (variables).
Each row (observation) is a unique tornado.
Look inside the simple feature data frame with the glimpse() function from the {dplyr} package.
dplyr::glimpse(Tornadoes.sf)## Rows: 63,645
## Columns: 23
## $ om <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ yr <dbl> 1950, 1950, 1950, 1950, 1950, 1950, 1950, 1950, 1950, 1950, 1…
## $ mo <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ dy <dbl> 3, 3, 3, 13, 25, 25, 26, 11, 11, 11, 11, 12, 12, 12, 12, 12, …
## $ date <chr> "1950-01-03", "1950-01-03", "1950-01-03", "1950-01-13", "1950…
## $ time <chr> "11:00:00", "11:55:00", "16:00:00", "05:25:00", "19:30:00", "…
## $ tz <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
## $ st <chr> "MO", "IL", "OH", "AR", "MO", "IL", "TX", "TX", "TX", "TX", "…
## $ stf <dbl> 29, 17, 39, 5, 29, 17, 48, 48, 48, 48, 48, 48, 48, 48, 48, 28…
## $ stn <dbl> 1, 2, 1, 1, 2, 3, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 10, 2, 1, …
## $ mag <dbl> 3, 3, 1, 3, 2, 2, 2, 2, 3, 2, 2, 2, 1, 2, 1, 2, 1, 3, 2, 4, 2…
## $ inj <dbl> 3, 3, 1, 1, 5, 0, 2, 0, 12, 5, 6, 8, 0, 0, 32, 2, 0, 15, 0, 7…
## $ fat <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 3, 0, 3, 0, 18, …
## $ loss <dbl> 6, 5, 4, 3, 5, 5, 0, 4, 4, 5, 5, 4, 4, 4, 5, 4, 0, 5, 3, 5, 5…
## $ closs <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ slat <dbl> 38.77, 39.10, 40.88, 34.40, 37.60, 41.17, 26.88, 29.42, 29.67…
## $ slon <dbl> -90.22, -89.30, -84.58, -94.37, -90.68, -87.33, -98.12, -95.2…
## $ elat <dbl> 38.8300, 39.1200, 40.8801, 34.4001, 37.6300, 41.1701, 26.8800…
## $ elon <dbl> -90.0300, -89.2300, -84.5799, -94.3699, -90.6500, -87.3299, -…
## $ len <dbl> 9.5, 3.6, 0.1, 0.6, 2.3, 0.1, 4.7, 9.9, 12.0, 4.6, 4.5, 8.0, …
## $ wid <dbl> 150, 130, 10, 17, 300, 100, 133, 400, 1000, 100, 67, 833, 233…
## $ fc <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ geometry <LINESTRING [°]> LINESTRING (-90.22 38.77, -..., LINESTRING (-89.3 …
The first 22 columns are variables (attributes). The last column contains the geometry. Information in the geometry column is in well-known text (WKT) format.
Each tornado is a coded as a LINESTRING with a start and end location. This is where the tm_shape() function looks for the geographic information.
Here you make a map showing the tracks of all the tornadoes since 2011. First filter the data frame keeping only tornadoes occurring after the year (yr) 2010.
TornadoesSince2011.sf <-
Tornadoes.sf |>
dplyr::filter(yr >= 2011) Next get a file containing the boundaries of the lower 48 states.
USA_48.sf <- USAboundaries::us_states() |>
dplyr::filter(!state_name %in% c("Hawaii", "Alaska", "Puerto Rico"))Then use the tm_shape() function together with the tm_borders() layer to draw the boundaries before adding the tornadoes. The tornadoes are in a separate spatial data frame so you use the tm_shape() function together with the tm_lines() layer.
tm_shape(shp = USA_48.sf, projection = 5070) +
tm_borders() +
tm_shape(shp = TornadoesSince2011.sf) +
tm_lines(col = "red")
The objects named TornadoesSince2011.sf and USA_48.sf are simple feature data frames. You map variables in the data frames as layers with successive calls to the tm_shape() function.
The default projection is geographic (latitude-longitude) which is changed using the projection = argument and specifying a EPSG number (or proj4 string). Here you use 5070 corresponding to USA Contiguous Albers Equal Area Conic, USGS (EPSG = 5070 or 102003).
You make the map interactive by first turning on the "view" mode with the tmap_mode() function before running the code.
tmap_mode("view")## tmap mode set to interactive viewing
tm_shape(USA_48.sf) +
tm_borders() +
tm_shape(TornadoesSince2011.sf) +
tm_lines(col = "red")